* Subroutine: M28Y12, M30Y12, M28BYS, M28CYS, M28CYN

      SUBROUTINE M28Y12(N, NZ, A, LICN, IRN, LIRN, ICN, U, IKEEP,IW, W,
     * LBLOCK, GROW, ABORT, IDISP, IFAIL, THR)
C     DERIVED FROM HARWELL LIBRARY ROUTINE MA28A
C     THROUGH F01BRF MARK 7 RELEASE. NAG 1978
C     AND Y12M -- ZLATEV & 0STERBY
C
C     THE PARAMETERS ARE AS FOLLOWS.....
C     N     INTEGER  ORDER OF MATRIX  NOT ALTERED BY SUBROUTINE.
C     NZ    INTEGER  NUMBER OF NON-ZEROS IN INPUT MATRIX  NOT
C     .     ALTERED BY SUBROUTINE.
C     A     REAL ARRAY  LENGTH LICN.  HOLDS NON-ZEROS OF
C     .     MATRIX ON ENTRY AND NON-ZEROS OF FACTORS ON EXIT.
C     .     REORDERED BY F01BRZ AND F01BRY AND ALTERED BY M30Y12.
C     LICN  INTEGER  LENGTH OF ARRAYS A AND ICN.  NOT ALTERED BY
C     .     SUBROUTINE.
C     IRN   INTEGER ARRAY  LENGTH LIRN.  HOLDS ROW INDICES ON INPUT
C     .     AND USED AS WORKSPACE BY M30Y12 TO HOLD COLUMN
C     .     ORIENTATION OF MATRIX.
C     LIRN  INTEGER  LENGTH OF ARRAY IRN.
C     ICN   INTEGER ARRAY  LENGTH LICN.  HOLDS COLUMN INDICES ON
C     .     ENTRY AND COLUMN INDICES OF DECOMPOSED MATRIX ON EXIT.
C     .     REORDERED BY F01BRZ AND F01BRY AND ALTERED BY M30Y12.
C     U     REAL VARIABLE  SET BY USER TO CONTROL BIAS TOWARDS
C     .     NUMERIC OR SPARSITY PIVOTING.  U=1.0 GIVES PARTIAL
C     .     PIVOTING WHILE U=0. DOES NOT CHECK MULTIPLIERS AT ALL.
C     .     VALUES OF U GREATER THAN ONE ARE TREATED AS ONE WHILE
C     .     NEGATIVE VALUES ARE TREATED AS ZERO.  NOT ALTERED BY
C     .     SUBROUTINE.
C     IKEEP  INTEGER ARRAY  LENGTH 5*N  USED AS WORKSPACE BY M28Y12
C     .     (SEE LATER COMMENTS).  IT IS NOT REQUIRED TO BE SET ON
C     .     ENTRY AND, ON EXIT, IT CONTAINS INFORMATION ABOUT THE
C     .     DECOMPOSITION.  IT SHOULD BE PRESERVED BETWEEN THIS CALL
C     .     AND SUBSEQUENT CALLS TO M28BSF OR F04AXF.
C     IKEEP(I,1),I=1,N  HOLDS THE TOTAL LENGTH OF THE PART OF ROW I
C     .     IN THE DIAGONAL BLOCK.
C     ROW IKEEP(I,2),I=1,N  OF THE INPUT MATRIX IS THE ITH ROW IN
C     .     PIVOT ORDER.
C     COLUMN IKEEP(I,3),I=1,N  OF THE INPUT MATRIX IS THE ITH
C     .     COLUMN IN PIVOT ORDER.
C     IKEEP(I,4),I=1,N  HOLDS THE LENGTH OF THE PART OF ROW I IN
C     .     THE L PART OF THE L/U DECOMPOSITION.
C     IKEEP(I,5),I=1,N  HOLDS THE LENGTH OF THE PART OF ROW I IN
C     .     THE OFF-DIAGONAL BLOCKS.  IF THERE IS ONLY ONE DIAGONAL
C     .     BLOCK, IKEEP(1,5) WILL BE SET TO -1.
C     IW    INTEGER ARRAY LENGTH 8*N.
C     .     USED AS WORKSPACE.
C     W     REAL ARRAY  LENGTH N.  USED BY F01BRQ BOTH AS WORKSPACE
C     .     AND TO RETURN GROWTH ESTIMATE IN W(1).  THE USE OF THIS
C     .     ARRAY BY M28Y12 IS THUS OPTIONAL DEPENDING ON THE VALUE
C     .     OF PARAMETER GROW.
C     LBLOCK  LOGICAL IF TRUE, F01BRY IS USED TO FIRST PERMUTE
C     .     THE MATRIX TO BLOCK LOWER TRIANGULAR FORM.
C     GROW  LOGICAL IF TRUE, THEN AN ESTIMATE OF THE INCREASE
C     .     IN SIZE OF MATRIX ELEMENTS DURING L/U DECOMPOSITION IS
C     .     GIVEN BY F01BRQ IN W(1).
C     ABORT  LOGICAL ARRAY LENGTH 4.
C     IF ABORT(1)=TRUE, THE ROUTINE WILL EXIT IMMEDIATELY
C     .     ON DETECTING STRUCTURAL SINGULARITY.
C     IF ABORT(2)=TRUE, THE ROUTINE WILL EXIT IMMEDIATELY
C     .     ON DETECTING NUMERICAL SINGULARITY.
C     IF ABORT(3)=TRUE, THE ROUTINE WILL EXIT IMMEDIATELY WHEN THE
C     .     AVAILABLE SPACE IN A/ICN IS FILLED UP BY THE PREVIOUSLY
C     .     DECOMPOSED ACTIVE, AND UNDECOMPOSED PARTS OF THE MATRIX.
C     IF ABORT(4)=TRUE, THE ROUTINE WILL EXIT IMMEDIATELY
C     .     ON DETECTING DUPLICATE ELEMENTS IN THE INPUT MATRIX.
C     IDISP  INTEGER ARRAY LENGTH 11. USED TO COMMUNICATE
C     .     INFORMATION ABOUT THE DECOMPOSITION TO THE USER AND
C     .     ALSO BETWEEN THIS CALL TO M28Y12 AND SUBSEQUENT CALLS
C     .     TO M28BSF AND F04AXF.
C     ON EXIT -
C     IDISP(1) AND IDISP(2) INDICATE THE POSITION IN ARRAYS A AND
C     .     ICN OF THE FIRST AND LAST ELEMENTS IN THE L/U
C     .     DECOMPOSITION OF THE DIAGONAL BLOCKS RESPECTIVELY.
C     IDISP(3)= IRNCP  THE NUMBER OF COMPRESSES ON ARRAY IRN.
C     IDISP(4)= ICNCP  THE NUMBER OF COMPRESSES ON ARRAYS ICN/A.
C     IDISP(5)= IRANK  ESTIMATED RANK OF THE MATRIX.
C     IDISP(6)= MINIRN MINIMUM LENGTH OF ARRAY IRN FOR SUCCESS ON
C     .     FUTURE RUNS.
C     IDISP(7)= MINICN MINIMUM LENGTH OF ARRAYS ICN/A FOR SUCCESS
C     .     ON FUTURE RUNS.
C     IDISP(8)= NUMNZ  STRUCTURAL RANK OF MATRIX.
C     IDISP(9)= NUM    NUMBER OF DIAGONAL BLOCKS.
C     IDISP(10)=LARGE  SIZE OF LARGEST DIAGONAL BLOCK.
C     IDISP(11)=Number of singular row.
C     IFAIL  INTEGER VARIABLE  USED AS ERROR FLAG BY ROUTINE.
C     THR    REAL VARIABLE  SET BY USER AS ZLATEV THRESHOLD
C
C     INTERNAL VARIABLES AND WORKSPACE USED IN  M28Y12 ARE DEFINED
C     WITHIN THE SUBROUTINE IMMEDIATELY PRIOR TO THEIR FIRST USE.
C     .. SCALAR ARGUMENTS ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='M28Y12')
      DOUBLE PRECISION U, THR
      INTEGER IFAIL, LICN, LIRN, N, NZ
      LOGICAL GROW, LBLOCK
C     .. ARRAY ARGUMENTS ..
      DOUBLE PRECISION A(LICN), W(N)
      INTEGER ICN(LICN), IDISP(10), IKEEP(N,5), IRN(LIRN), IW(N,8)
      LOGICAL ABORT(4)
C     ..
C     .. LOCAL SCALARS ..
C$P 1
      DOUBLE PRECISION THEMAX, UPRIV, ZERO
      INTEGER I1, I, IEND, II, IPRIV4, ISAVE, J1, J2, J, JAY, JJ,KNUM,
     * LENGTH, MOVE, NADV, NERR, NEWJ1, NEWPOS
C     .. Local Arrays ..
      CHARACTER*1       P01REC(1)
      CHARACTER*65      REC(3)
C     .. External Functions ..
      INTEGER           P01ABF
      EXTERNAL          P01ABF
C     .. External Subroutines ..
      EXTERNAL          F01BRQ, F01BRR, M30Y12, F01BRY, F01BRZ, X04AAF,
     *                  X04ABF, X04BAF
C     .. Intrinsic Functions ..
c      INTRINSIC         ABS, MAX, MOD !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C     .. Data statements ..
      DATA ZERO /0.0D0/
      ISAVE = IFAIL
      IFAIL = 0
C     NERR IS THE UNIT NUMBER FOR ERROR MESSAGES
C     NADV IS THE UNIT NUMBER FOR DUPLICATE ELEMENT WARNING MESSAGES
      CALL X04AAF(0, NERR)
      CALL X04ABF(0, NADV)
C     UPRIV PRIVATE COPY OF U IS USED IN CASE IT IS OUTSIDE
C     RANGE  ZERO TO ONE  AND  IS THUS ALTERED BY M30Y12.
      UPRIV = U
C     SIMPLE DATA CHECK ON INPUT VARIABLES AND ARRAY DIMENSIONS.
      IF (N.GT.0) GO TO 20
      IFAIL = 8
      GO TO 320
   20 IF (NZ.GT.0) GO TO 40
      IFAIL = 9
      GO TO 320
   40 IF (LICN.GE.NZ) GO TO 60
      IFAIL = 10
      GO TO 320
   60 IF (LIRN.GE.NZ) GO TO 80
      IFAIL = 11
      GO TO 320
C
C     DATA CHECK TO SEE IF ALL INDICES LIE BETWEEN 1 AND N.
   80 DO 100 I=1,NZ
         IF (IRN(I).GT.0 .AND. IRN(I).LE.N .AND. ICN(I).GT.0 .AND.ICN(I)
     *   .LE.N) GO TO 100
C     ** CODE FOR OUTPUT OF ERROR MESSAGE **************************
         IF (MOD(ISAVE/10,10).NE.0) THEN
            WRITE (REC,FMT=99999) I, A(I), IRN(I), ICN(I)
            CALL X04BAF(NERR,REC(1))
            CALL X04BAF(NERR,REC(2))
            CALL X04BAF(NERR,REC(3))
         END IF
C     ** END OF CODE FOR OUTPUT OF ERROR MESSAGE *******************
         IFAIL = 12
  100 CONTINUE
      IF (IFAIL.GT.0) GO TO 340
C
C     SORT ELEMENTS INTO ROW ORDER.
C
      CALL F01BRZ(N, NZ, A, ICN, IW, IRN)
      IPRIV4 = IW(1,1)
C     PART OF IKEEP IS USED HERE AS A WORK-ARRAY.  IKEEP(I,2) IS
C     THE LAST ROW TO HAVE A NON-ZERO IN COLUMN I.  IKEEP(I,3)
C     IS THE OFF-SET OF COLUMN I FROM THE START OF THE ROW.
      DO 120 I=1,N
         IKEEP(I,2) = 0
         IKEEP(I,1) = 0
  120 CONTINUE
C
C     CHECK FOR DUPLICATE ELEMENTS .. SUMMING ANY SUCH ENTRIES AND
C     PRINTING A WARNING MESSAGE ON UNIT NADV.
C     MOVE IS EQUAL TO THE NUMBER OF DUPLICATE ELEMENTS FOUND.
      MOVE = 0
C     THE LOOP ALSO CALCULATES THE LARGEST ELEMENT IN THE MATRIX,
C     THEMAX.
      THEMAX = ZERO
C     J1 IS POSITION IN ARRAYS OF FIRST NON-ZERO IN ROW.
      J1 = IPRIV4
      DO 220 I=1,N
         IF (I.NE.N) GO TO 140
         IPRIV4 = NZ + 1
         GO TO 160
  140    IPRIV4 = IW(I+1,1)
  160    LENGTH = IPRIV4 - J1
         IF (LENGTH.EQ.0) GO TO 220
         J2 = IPRIV4 - 1
         NEWJ1 = J1 - MOVE
         DO 200 JJ=J1,J2
            J = ICN(JJ)
            THEMAX = DMAX1(THEMAX,DABS(A(JJ)))
            IF (IKEEP(J,2).EQ.I) GO TO 180
C     FIRST TIME COLUMN HAS OCURRED IN CURRENT ROW.
            IKEEP(J,2) = I
            IKEEP(J,3) = JJ - MOVE - NEWJ1
            IF (MOVE.EQ.0) GO TO 200
C     SHIFT NECESSARY BECAUSE OF  PREVIOUS DUPLICATE ELEMENT.
            NEWPOS = JJ - MOVE
            A(NEWPOS) = A(JJ)
            ICN(NEWPOS) = ICN(JJ)
            GO TO 200
C     DUPLICATE ELEMENT.
  180       MOVE = MOVE + 1
            LENGTH = LENGTH - 1
            JAY = IKEEP(J,3) + NEWJ1
C     ** CODE FOR OUTPUT OF WARNING MESSAGE ************************
            IF (MOD(ISAVE/100,100).NE.0) THEN
               WRITE (REC,FMT=99998) I, J, A(JJ)
               CALL X04BAF(NADV,REC(1))
               CALL X04BAF(NADV,REC(2))
               CALL X04BAF(NADV,REC(3))
            END IF
C     ** END OF CODE FOR OUTPUT OF WARNING MESSAGE *****************
            IF (ABORT(4)) IFAIL = 13
            A(JAY) = A(JAY) + A(JJ)
            THEMAX = DMAX1(THEMAX,DABS(A(JAY)))
  200    CONTINUE
         IKEEP(I,1) = LENGTH
         J1 = IPRIV4
  220 CONTINUE
C     IF ABORT(4) IS .TRUE., DUPLICATE ELEMENTS ARE REGARDED AS
C     AN ERROR
      IF (IFAIL.GT.0) GO TO 320
C
C     KNUM IS ACTUAL NUMBER OF NON-ZEROS IN MATRIX WITH ANY
C     MULTIPLE ENTRIES COUNTED ONLY ONCE.
      KNUM = NZ - MOVE
      IF (.NOT.LBLOCK) GO TO 240
C
C     PERFORM BLOCK TRIANGULARISATION.
      IFAIL = ISAVE
      CALL F01BRY(N, ICN, A, LICN, IKEEP, IDISP, IKEEP(1,2),IKEEP(1,3),
     * IKEEP(1,5), IW(1,3), IW, ABORT(1), IFAIL)
C     ON EXIT FROM F01BRY IDISP(8)=NUMNZ, IDISP(9)=NUM,
C     IDISP(10) = LARGE.
      IF (IFAIL.NE.0) GO TO 340
      GO TO 300
C
C     BLOCK TRIANGULARIZATION NOT REQUESTED.
C     MOVE STRUCTURE TO END OF DATA ARRAYS IN PREPARATION FOR
C     M30Y12.
C     ALSO SET LENOFF(1) TO -1 AND SET PERMUTATION ARRAYS.
  240 DO 260 I=1,KNUM
         II = KNUM - I + 1
         NEWPOS = LICN - I + 1
         ICN(NEWPOS) = ICN(II)
         A(NEWPOS) = A(II)
  260 CONTINUE
      IDISP(1) = 1
      IDISP(2) = LICN - KNUM + 1
      DO 280 I=1,N
         IKEEP(I,2) = I
         IKEEP(I,3) = I
  280 CONTINUE
      IKEEP(1,5) = -1
C
C     PERFORM L/U DECOMOSITION ON DIAGONAL BLOCKS.
  300 IFAIL = ISAVE
      CALL M30Y12(N, ICN, A, LICN, IKEEP, IKEEP(1,4), IDISP,IKEEP(1,2),
     * IKEEP(1,3), IRN, LIRN, IW(1,2), IW(1,3),IW(1,4), IW(1,5), IW(1,6)
     *, IW(1,7), IW(1,8), IW, UPRIV,ABORT, IFAIL, THR)
C
C     ON EXIT FROM M30Y12 IDISP(3) TO IDISP(7) CONTAIN
C     IRNCP, ICNCP, IRANK, MINIRN AND MINICN RESPECTIVELY.
      IDISP(6) = MAX0(IDISP(6),NZ)
      IDISP(7) = MAX0(IDISP(7),NZ)
      IF (IFAIL.GT.0) GO TO 340
C
C     REORDER OFF-DIAGONAL BLOCKS ACCORDING TO PIVOT PERMUTATION.
      I1 = IDISP(1) - 1
      IF (I1.NE.0) CALL F01BRR(N, ICN, A, I1, IKEEP(1,5),IKEEP(1,2),
     * IKEEP(1,3), IW, IRN)
C
C     OPTIONALLY CALCULATE ELEMENT GROWTH ESTIMATE.
      I1 = IDISP(1)
      IEND = LICN - I1 + 1
      IF (GROW) CALL F01BRQ(N, ICN, A(I1), IEND, IKEEP, IKEEP(1,4),W)
C     INCREMENT GROWTH ESTIMATE BY ORIGINAL MAXIMUM ELEMENT.
      IF (GROW) W(1) = W(1) + THEMAX
      IF (IFAIL.EQ.0) GO TO 360
C
  320 CONTINUE
C     ** CODE FOR OUTPUT OF ERROR MESSAGES *************************
      IF (MOD(ISAVE/10,10).EQ.0) GO TO 340
      IF (IFAIL.EQ.-2) WRITE (REC,FMT=99992)
      IF (IFAIL.EQ.-1) WRITE (REC,FMT=99991)
      IF (IFAIL.EQ.8) WRITE (REC,FMT=99997) N
      IF (IFAIL.EQ.9) WRITE (REC,FMT=99996) NZ
      IF (IFAIL.EQ.10) WRITE (REC,FMT=99995) LICN, NZ
      IF (IFAIL.EQ.11) WRITE (REC,FMT=99994) LIRN, NZ
      IF (IFAIL.EQ.13) WRITE (REC,FMT=99993)
      IF (IFAIL.EQ.-1 .OR. IFAIL.EQ.-2 .OR. IFAIL.EQ.13 .OR.
     *    (IFAIL.GE.8 .AND. IFAIL.LE.11)) THEN
         CALL X04BAF(NERR,REC(1))
         CALL X04BAF(NERR,REC(2))
      END IF
C     ** END OF CODE FOR OUTPUT OF ERROR MESSAGES ******************
C
  340 IFAIL = P01ABF(ISAVE,IFAIL,SRNAME,0,P01REC)
  360 RETURN
99999 FORMAT (/' ON ENTRY IRN(I) OR ICN(I) IS OUT OF RANGE - I = ',I8,
     *  /'  A(I) = ',1P,D14.6,'  IRN(I) = ',I8,'  ICN(I) = ',I8)
99998 FORMAT (/' M28Y12 FOUND DUPLICATE ELEMENT WITH INDICES ',I6,', ',
     *  I6,/'  VALUE = ',1P,D14.6)
99997 FORMAT (/' ON ENTRY N .LE. 0 , N = ',I10)
99996 FORMAT (/' ON ENTRY NZ .LE. 0 , NZ = ',I10)
99995 FORMAT (/' ON ENTRY LICN .LT. NZ , LICN = ',I8,'  NZ = ',I8)
99994 FORMAT (/' ON ENTRY LIRN .LT. NZ , LIRN = ',I8,'  NZ = ',I8)
99993 FORMAT (/' DUPLICATE ELEMENTS FOUND ON INPUT - SEE ADVISORY MESS',
     *  'AGES')
99992 FORMAT (/' MATRIX IS NUMERICALLY SINGULAR - DECOMPOSITION COMPLE',
     *  'TED')
99991 FORMAT (/' MATRIX IS STRUCTURALLY SINGULAR - DECOMPOSITION COMPL',
     *  'ETED')
      END

      SUBROUTINEM30Y12(NN,ICN,A,LICN,LENR,LENRL,IDISP,IP,IQ,IRN,LIRN,LEN 
     *C,IFIRST,LASTR,NEXTR,LASTC,NEXTC,IPTR,IPC,U,ABORT,IFAIL,THR) 
      DOUBLEPRECISIONU 
      DOUBLEPRECISIONTHR 
      INTEGERIFAIL,LICN,LIRN,NN 
      DOUBLEPRECISIONA(LICN) 
      INTEGERICN(LICN),IDISP(11),IFIRST(NN),IP(NN),IPC(NN),IPTR(NN),IQ(N 
     *N),IRN(LIRN),LASTC(NN),LASTR(NN),LENC(NN),LENR(NN),LENRL(NN),NEXTC 
     *(NN),NEXTR(NN) 
      LOGICALABORT(3) 
      DOUBLEPRECISIONRMAX,AU,UMAX,ZERO 
      INTEGERDISPC,I1,I2,I,IACTIV,IBEG,ICNCP,IDISPC,IDUMMY,IEND,IFILL,IF 
     *IR,II,III,IJFIR,IJP1,IJPOS,ILAST,INDROW,IOP,IPIV,IPOS,IRANK,IRNCP, 
     *IROWS,ISAVE,ISING,ISTART,ISW1,ISW,ITOP,J1,J2,J,JBEG,JCOST,JCOUNT,J 
     *DIFF,JEND,JJ,JMORE,JNPOS,JOLD,JPIV,JPOS,JROOM,JVAL,JZER,JZERO,K,KC 
     *OST,L,LC,LENPIV,LL,LR,MINICN,MINIRN,MOREI,N,NBLOCK,NC,NERR,NNM1,NR 
     *,NUM,NZ2,NZ,NZCOL,NZMIN,NZPC,NZROW,OLDEND,OLDPIV,PIVEND,PIVOT,PIVR 
     *OW,ROWI 
      DATANPZLAT/3/ 
      DATAUMAX/.9999D0/ 
      DATAZERO/0.0D0/ 
      ISAVE=IFAIL 
      IFAIL=0 
      CALLX04AAF(0,NERR) 
      MINIRN=0 
      MINICN=IDISP(1)-1 
      MOREI=0 
      IRANK=NN 
      IRNCP=0 
      ICNCP=0 
      U=DMIN1(U,UMAX) 
      U=DMAX1(U,ZERO) 
      IBEG=IDISP(1) 
      IACTIV=IDISP(2) 
      NZROW=LICN-IACTIV+1 
      MINICN=NZROW+MINICN 
      NUM=1 
      IPTR(1)=IACTIV 
      IF(.NOT.((NN.NE.1)))GOTO32757 
      NNM1=NN-1 
      DO32754I=1,NNM1 
      IF(IP(I).LT.0)NUM=NUM+1 
      IPTR(I+1)=IPTR(I)+LENR(I) 
32754 CONTINUE 
32757 CONTINUE 
      ILAST=0 
      DO32751NBLOCK=1,NUM 
      ISTART=ILAST+1 
      DO32748IROWS=ISTART,NN 
      IF(IP(IROWS).LT.0)GOTO80 
32748 CONTINUE 
      IROWS=NN 
80    ILAST=IROWS 
      N=ILAST-ISTART+1 
      IF(.NOT.(N.EQ.1))GOTO32745 
      LENRL(ILAST)=0 
      ISING=ISTART 
      IF(.NOT.((LENR(ILAST).EQ.0)))GOTO32742 
      IRANK=IRANK-1 
      ISING=-ISING 
      IF(IFAIL.NE.-2.AND.IFAIL.NE.5)IFAIL=-1 
      IF(.NOT.((ABORT(1))))GOTO32739 
      IFAIL=1 
      GOTO2020 
32739 CONTINUE 
      GOTO32741 
32742 CONTINUE 
      IF(.NOT.((A(IACTIV).EQ.ZERO)))GOTO32736 
      ISING=-ISING 
      IRANK=IRANK-1 
      IPTR(ILAST)=0 
      IF(IFAIL.NE.5)IFAIL=-2 
      IF(.NOT.((ABORT(2))))GOTO32733 
      IFAIL=2 
      GOTO2020 
32733 CONTINUE 
      GOTO32735 
32736 CONTINUE 
      A(IBEG)=A(IACTIV) 
      ICN(IBEG)=ICN(IACTIV) 
      IACTIV=IACTIV+1 
      IPTR(ISTART)=0 
      IBEG=IBEG+1 
      NZROW=NZROW-1 
32735 CONTINUE 
32741 CONTINUE 
      LASTR(ISTART)=ISTART 
      LASTC(ISTART)=ISING 
      GOTO32744 
32745 CONTINUE 
      IF(.NOT.(ILAST.EQ.NN))GOTO32730 
      ITOP=LICN 
      GOTO32729 
32730 CONTINUE 
      ITOP=IPTR(ILAST+1)-1 
32729 CONTINUE 
      DO32727I=ISTART,ILAST 
      LENRL(I)=0 
      LENC(I)=0 
32727 CONTINUE 
      IF(.NOT.((ITOP-IACTIV.GE.LIRN)))GOTO32724 
      MINIRN=ITOP-IACTIV+1 
      PIVOT=ISTART-1 
      GOTO1980 
32724 CONTINUE 
      DO32721II=IACTIV,ITOP 
      I=ICN(II) 
      LENC(I)=LENC(I)+1 
32721 CONTINUE 
      IPC(ILAST)=LIRN+1 
      J1=ISTART+1 
      DO32718JJ=J1,ILAST 
      J=ILAST-JJ+J1-1 
      IPC(J)=IPC(J+1)-LENC(J+1) 
32718 CONTINUE 
      DO32715INDROW=ISTART,ILAST 
      J1=IPTR(INDROW) 
      J2=J1+LENR(INDROW)-1 
      IF(.NOT.((J1.LE.J2)))GOTO32712 
      DO32709JJ=J1,J2 
      J=ICN(JJ) 
      IPOS=IPC(J)-1 
      IRN(IPOS)=INDROW 
      IPC(J)=IPOS 
32709 CONTINUE 
32712 CONTINUE 
32715 CONTINUE 
      DISPC=IPC(ISTART) 
      NZCOL=LIRN-DISPC+1 
      MINIRN=MAX0(NZCOL,MINIRN) 
      NZMIN=1 
      DO32706INZ=1,N 
      IFIRST(INZ)=0 
32706 CONTINUE 
      DO32703II=ISTART,ILAST 
      I=ILAST-II+ISTART 
      NZ=LENR(I) 
      IF(.NOT.(NZ.EQ.0))GOTO32700 
      IPTR(I)=0 
      LASTR(I)=0 
      GOTO32699 
32700 CONTINUE 
      ISW=IFIRST(NZ) 
      IFIRST(NZ)=I 
      IF(.NOT.(ISW.LE.0))GOTO32697 
      NEXTR(I)=0 
      LASTR(I)=ISW 
      GOTO32696 
32697 CONTINUE 
      NEXTR(I)=ISW 
      LASTR(I)=LASTR(ISW) 
      LASTR(ISW)=I 
32696 CONTINUE 
32699 CONTINUE 
32703 CONTINUE 
      DO32694PIVOT=ISTART,ILAST 
      NSCAN=MIN0(NPZLAT,ILAST+1-PIVOT) 
      NZ2=NZMIN 
      JCOST=N*N 
      DO660L=1,2 
      LL=L 
      ISCAN=0 
      NZ=NZ2 
32691 IF(.NOT.(ISCAN.LT.NSCAN.AND.NZ.LE.N))GOTO32690 
      IJFIR=IFIRST(NZ) 
      IF(.NOT.(IJFIR.GT.0))GOTO32688 
      LL=2 
32685 IF(.NOT.(ISCAN.LT.NSCAN.AND.IJFIR.NE.0))GOTO32684 
      ISCAN=ISCAN+1 
      I=IJFIR 
      IJFIR=NEXTR(I) 
      RMAX=ZERO 
      J1=IPTR(I)+LENRL(I) 
      IF(.NOT.(ISCAN.EQ.1))GOTO32682 
      IJPOS=J1 
      IPIV=I 
      JPIV=ICN(J1) 
32682 CONTINUE 
      J2=IPTR(I)+LENR(I)-1 
      DO32679JJ=J1,J2 
      RMAX=DMAX1(RMAX,DABS(A(JJ))) 
32679 CONTINUE 
      AU=RMAX*U 
      DO32676JJ=J1,J2 
      IF(.NOT.((DABS(A(JJ)).GT.AU.OR.L.NE.1)))GOTO32673 
      J=ICN(JJ) 
      KCOST=(NZ-1)*(LENC(J)-1) 
      IF(.NOT.((KCOST.LT.JCOST).OR.(KCOST.EQ.JCOST.AND.DABS(A(JJ)).GT.DA 
     *BS(A(IJPOS)))))GOTO32670 
      JCOST=KCOST 
      IJPOS=JJ 
      IPIV=I 
      JPIV=J 
32670 CONTINUE 
32673 CONTINUE 
32676 CONTINUE 
      GOTO32685 
32684 CONTINUE 
      GOTO32687 
32688 CONTINUE 
      IF(.NOT.(LL.EQ.1.AND.IJFIR.EQ.0))GOTO32667 
      NZMIN=NZ+1 
32667 CONTINUE 
32687 CONTINUE 
      NZ=NZ+1 
      GOTO32691 
32690 CONTINUE 
      IF(ISCAN.EQ.NSCAN)GOTO820 
      IRANK=IRANK-1 
660   CONTINUE 
      WRITE(5,662)ISCAN 
662   Format(' ISCAN=',I8) 
      IF(IFAIL.NE.-2.AND.IFAIL.NE.5)IFAIL=-1 
      IRANK=IRANK-ILAST+PIVOT+1 
      IF(.NOT.(ABORT(1)))GOTO32664 
      IFAIL=1 
      GOTO2020 
32664 CONTINUE 
      K=PIVOT-1 
      I=ISTART 
32661 IF(.NOT.(I.LE.ILAST.AND.LENRL(I).NE.0))GOTO32660 
      IF(.NOT.(LASTR(I).EQ.0))GOTO32658 
      K=K+1 
      LASTR(I)=K 
      IF(.NOT.(LENRL(I).NE.0))GOTO32655 
      MINICN=MAX0(MINICN,NZROW+IBEG-1+MOREI+LENRL(I)) 
      IF(.NOT.(IACTIV-IBEG.LT.LENRL(I)))GOTO32652 
      CALLF01BRS(A,ICN,IPTR(ISTART),N,IACTIV,ITOP,.TRUE.,ICNCP) 
      IF(.NOT.(IACTIV-IBEG.LT.LENRL(I)))GOTO32649 
      MOREI=MOREI+IBEG-IDISP(1) 
      IBEG=IDISP(1) 
      IFAIL=5 
      IF(ABORT(3))GOTO2020 
      IF(MOD(ISAVE/10,10).NE.0)WRITE(NERR,99997) 
32649 CONTINUE 
32652 CONTINUE 
      J1=IPTR(I) 
      J2=J1+LENRL(I)-1 
      IPTR(I)=0 
      DO32646JJ=J1,J2 
      A(IBEG)=A(JJ) 
      ICN(IBEG)=ICN(JJ) 
      ICN(JJ)=0 
      IBEG=IBEG+1 
32646 CONTINUE 
      NZROW=NZROW-LENRL(I) 
32655 CONTINUE 
32658 CONTINUE 
      I=I+1 
      GOTO32661 
32660 CONTINUE 
780   K=PIVOT-1 
      DO32643I=ISTART,ILAST 
      IF(.NOT.((LASTC(I).EQ.0)))GOTO32640 
      K=K+1 
      LASTC(I)=-K 
      IF(K.EQ.ILAST)GOTO32693 
32640 CONTINUE 
32643 CONTINUE 
820   ISING=PIVOT 
      IF(.NOT.((A(IJPOS).EQ.ZERO)))GOTO32637 
      Idisp(11)=IP(IPIV) 
      Write(6,*)' I sing=',IP(IPIV) 
      ISING=-ISING 
      IF(IFAIL.NE.5)IFAIL=-2 
      IF(.NOT.((ABORT(2))))GOTO32634 
      IFAIL=2 
      GOTO2020 
32634 CONTINUE 
32637 CONTINUE 
      OLDPIV=IPTR(IPIV)+LENRL(IPIV) 
      OLDEND=IPTR(IPIV)+LENR(IPIV)-1 
      I1=IPC(JPIV) 
      I2=I1+LENC(JPIV)-1 
      DO32631II=I1,I2 
      I=IRN(II) 
      LR=LASTR(I) 
      NR=NEXTR(I) 
      IF(NR.NE.0)LASTR(NR)=LR 
      IF(.NOT.(LR.GT.0))GOTO32628 
      NEXTR(LR)=NR 
      GOTO32627 
32628 CONTINUE 
      NZ=LENR(I)-LENRL(I) 
      IF(.NOT.(NR.NE.0))GOTO32625 
      IFIRST(NZ)=NR 
      GOTO32624 
32625 CONTINUE 
      IFIRST(NZ)=0 
32624 CONTINUE 
32627 CONTINUE 
32631 CONTINUE 
      LASTC(JPIV)=ISING 
      LASTR(IPIV)=PIVOT 
      IF(.NOT.(OLDPIV.NE.IJPOS))GOTO32622 
      AU=A(OLDPIV) 
      A(OLDPIV)=A(IJPOS) 
      A(IJPOS)=AU 
      ICN(IJPOS)=ICN(OLDPIV) 
      ICN(OLDPIV)=JPIV 
32622 CONTINUE 
      MINICN=MAX0(MINICN,NZROW+IBEG-1+MOREI+LENR(IPIV)) 
      IF(.NOT.(IACTIV-IBEG.LT.LENR(IPIV)))GOTO32619 
      CALLF01BRS(A,ICN,IPTR(ISTART),N,IACTIV,ITOP,.TRUE.,ICNCP) 
      OLDPIV=IPTR(IPIV)+LENRL(IPIV) 
      OLDEND=IPTR(IPIV)+LENR(IPIV)-1 
      IF(.NOT.(IACTIV-IBEG.LT.LENR(IPIV)))GOTO32616 
      MOREI=MOREI+IBEG-IDISP(1) 
      IBEG=IDISP(1) 
      IFAIL=5 
      IF(ABORT(3))GOTO2020 
      IF(MOD(ISAVE/10,10).NE.0)WRITE(NERR,99997) 
      IF(.NOT.(IACTIV-IBEG.LT.LENR(IPIV)))GOTO32613 
      IFAIL=4 
      GOTO2020 
32613 CONTINUE 
32616 CONTINUE 
32619 CONTINUE 
      IJPOS=0 
      J1=IPTR(IPIV) 
      DO32610JJ=J1,OLDEND 
      A(IBEG)=A(JJ) 
      ICN(IBEG)=ICN(JJ) 
      IF(.NOT.((IJPOS.EQ.0)))GOTO32607 
      IF(ICN(JJ).EQ.JPIV)IJPOS=IBEG 
      ICN(JJ)=0 
      GOTO32606 
32607 CONTINUE 
      K=IBEG-IJPOS 
      J=ICN(JJ) 
      ICN(JJ)=IQ(J) 
      IQ(J)=-K 
32606 CONTINUE 
      IBEG=IBEG+1 
32610 CONTINUE 
      IJP1=IJPOS+1 
      PIVEND=IBEG-1 
      LENPIV=PIVEND-IJPOS 
      NZROW=NZROW-LENRL(IPIV)-1 
      IPTR(IPIV)=OLDPIV+1 
      IF(LENPIV.EQ.0)IPTR(IPIV)=0 
      DO32604JJ=IJPOS,PIVEND 
      J=ICN(JJ) 
      I1=IPC(J) 
      LENC(J)=LENC(J)-1 
      I2=IPC(J)+LENC(J)-1 
      IF(.NOT.((I2.GE.I1)))GOTO32601 
      DO32598II=I1,I2 
      IF(.NOT.((IRN(II).EQ.IPIV)))GOTO32595 
      IRN(II)=IRN(I2+1) 
      GOTO32597 
32595 CONTINUE 
32598 CONTINUE 
32597 CONTINUE 
32601 CONTINUE 
      IRN(I2+1)=0 
32604 CONTINUE 
      NZCOL=NZCOL-LENPIV-1 
      NZPC=LENC(JPIV) 
      IF(.NOT.(NZPC.NE.0))GOTO32592 
      DO32589III=1,NZPC 
      II=IPC(JPIV)+III-1 
      I=IRN(II) 
      J1=IPTR(I)+LENRL(I) 
      IEND=IPTR(I)+LENR(I)-1 
      JJ=J1 
32586 IF(.NOT.(ICN(JJ).NE.JPIV.AND.JJ.LE.IEND))GOTO32585 
      JJ=JJ+1 
      GOTO32586 
32585 CONTINUE 
      IF(.NOT.(JJ.LE.IEND))GOTO32583 
      IF(.NOT.(A(IJPOS).EQ.ZERO))GOTO32580 
      AU=ZERO 
      GOTO32579 
32580 CONTINUE 
      AU=-A(JJ)/A(IJPOS) 
32579 CONTINUE 
      A(JJ)=A(J1) 
      A(J1)=AU 
      ICN(JJ)=ICN(J1) 
      ICN(J1)=JPIV 
      LENRL(I)=LENRL(I)+1 
32583 CONTINUE 
      IF(.NOT.(LENPIV.NE.0))GOTO32577 
      ROWI=J1+1 
      IOP=0 
      IF(.NOT.(ROWI.LE.IEND))GOTO32574 
      DO1140JJ=ROWI,IEND 
      J=ICN(JJ) 
      IF(IQ(J).GT.0)GOTO1140 
      IOP=IOP+1 
      PIVROW=IJPOS-IQ(J) 
      A(JJ)=A(JJ)+AU*A(PIVROW) 
      ICN(PIVROW)=-ICN(PIVROW) 
1140  CONTINUE 
32574 CONTINUE 
      IFILL=LENPIV-IOP 
      DO1162JJ=IJP1,PIVEND 
      J=ICN(JJ) 
      IF(J.LT.0)GOTO1162 
      IF(.NOT.(DABS(AU*A(JJ)).LT.THR))GOTO32571 
      ICN(JJ)=-J 
      IFILL=IFILL-1 
32571 CONTINUE 
1162  CONTINUE 
      IF(.NOT.(IFILL.NE.0))GOTO32568 
      MINICN=MAX0(MINICN,MOREI+IBEG-1+NZROW+IFILL+LENR(I)) 
      DO32565JDIFF=1,IFILL 
      JNPOS=IEND+JDIFF 
      IF(JNPOS.GT.LICN)GOTO1200 
      IF(ICN(JNPOS).NE.0)GOTO1200 
32565 CONTINUE 
      IEND=IEND+1 
      GOTO1360 
1200  JMORE=IFILL-JDIFF+1 
      I1=IPTR(I) 
      DO32562JDIFF=1,JMORE 
      JNPOS=I1-JDIFF 
      IF(JNPOS.LT.IACTIV)GOTO1240 
      IF(ICN(JNPOS).NE.0)GOTO1260 
32562 CONTINUE 
1240  JNPOS=I1-JMORE 
      GOTO1280 
1260  JNPOS=IACTIV-LENR(I)-IFILL 
1280  continue 
      IF(.NOT.(JNPOS.LT.IBEG))GOTO32559 
      CALLF01BRS(A,ICN,IPTR(ISTART),N,IACTIV,ITOP,.TRUE.,ICNCP) 
      I1=IPTR(I) 
      IEND=I1+LENR(I)-1 
      JNPOS=IACTIV-LENR(I)-IFILL 
      IF(.NOT.(JNPOS.LT.IBEG))GOTO32556 
      MOREI=MOREI+IBEG-IDISP(1)-LENPIV-1 
      IFAIL=5 
      IF(ABORT(3))GOTO2020 
      IF(MOD(ISAVE/10,10).NE.0)WRITE(NERR,99997) 
      IBEG=IDISP(1) 
      ICN(IBEG)=JPIV 
      A(IBEG)=A(IJPOS) 
      IJPOS=IBEG 
      DO32553JJ=IJP1,PIVEND 
      IBEG=IBEG+1 
      A(IBEG)=A(JJ) 
      ICN(IBEG)=ICN(JJ) 
32553 CONTINUE 
      IJP1=IJPOS+1 
      PIVEND=IBEG 
      IBEG=IBEG+1 
      IF(.NOT.(JNPOS.LT.IBEG))GOTO32550 
      IFAIL=4 
      GOTO2020 
32550 CONTINUE 
32556 CONTINUE 
32559 CONTINUE 
      IACTIV=MIN0(IACTIV,JNPOS) 
      IPTR(I)=JNPOS 
      DO32547JJ=I1,IEND 
      A(JNPOS)=A(JJ) 
      ICN(JNPOS)=ICN(JJ) 
      JNPOS=JNPOS+1 
      ICN(JJ)=0 
32547 CONTINUE 
      IEND=JNPOS 
32568 CONTINUE 
1360  NZROW=NZROW+IFILL 
      IENBEG=IEND 
      DO32544JJ=IJP1,PIVEND 
      J=ICN(JJ) 
      IF(.NOT.(J.GT.0))GOTO32541 
      A(IEND)=AU*A(JJ) 
      ICN(IEND)=J 
      IEND=IEND+1 
      MINIRN=MAX0(MINIRN,NZCOL+LENC(J)+1) 
      JEND=IPC(J)+LENC(J) 
      JROOM=NZPC-III+1+LENC(J) 
      IF(.NOT.((JEND.LE.LIRN.AND.IRN(JEND).NE.0).OR.JEND.GT.LIRN))GOTO32 
     *538 
      IF(.NOT.(JROOM.GE.DISPC))GOTO32535 
      CALLF01BRS(A,IRN,IPC(ISTART),N,DISPC,LIRN,.FALSE.,IRNCP) 
      IF(.NOT.(JROOM.GE.DISPC))GOTO32532 
      JROOM=DISPC-1 
      IF(.NOT.(JROOM.LT.LENC(J)+1))GOTO32529 
      GOTO1980 
32529 CONTINUE 
32532 CONTINUE 
32535 CONTINUE 
      JBEG=IPC(J) 
      JEND=IPC(J)+LENC(J)-1 
      JZERO=DISPC-1 
      DISPC=DISPC-JROOM 
      IDISPC=DISPC 
      DO32526II=JBEG,JEND 
      IRN(IDISPC)=IRN(II) 
      IRN(II)=0 
      IDISPC=IDISPC+1 
32526 CONTINUE 
      IPC(J)=DISPC 
      JEND=IDISPC 
      DO32523II=JEND,JZERO 
      IRN(II)=0 
32523 CONTINUE 
32538 CONTINUE 
      IRN(JEND)=I 
      NZCOL=NZCOL+1 
      LENC(J)=LENC(J)+1 
      GOTO32540 
32541 CONTINUE 
      ICN(JJ)=-J 
32540 CONTINUE 
32544 CONTINUE 
      LENR(I)=LENR(I)+IFILL 
32577 CONTINUE 
32589 CONTINUE 
      I1=IPC(JPIV) 
      I2=IPC(JPIV)+LENC(JPIV)-1 
      NZCOL=NZCOL-LENC(JPIV) 
      DO32520II=I1,I2 
      I=IRN(II) 
      IRN(II)=0 
      NZ=LENR(I)-LENRL(I) 
      IF(.NOT.(NZ.EQ.0))GOTO32517 
      LASTR(I)=0 
      GOTO32516 
32517 CONTINUE 
      IFIR=IFIRST(NZ) 
      IFIRST(NZ)=I 
      IF(.NOT.(IFIR.NE.0))GOTO32514 
      LASTR(I)=LASTR(IFIR) 
      NEXTR(I)=IFIR 
      LASTR(IFIR)=I 
      GOTO32513 
32514 CONTINUE 
      LASTR(I)=0 
      NEXTR(I)=0 
      NZMIN=MIN0(NZMIN,NZ) 
32513 CONTINUE 
32516 CONTINUE 
32520 CONTINUE 
32592 CONTINUE 
      IPC(JPIV)=0 
      IF(.NOT.((LENPIV.NE.0)))GOTO32511 
      NZROW=NZROW-LENPIV 
      JVAL=IJP1 
      JZER=IPTR(IPIV) 
      IPTR(IPIV)=0 
      DO32508JCOUNT=1,LENPIV 
      J=ICN(JVAL) 
      IQ(J)=ICN(JZER) 
      ICN(JZER)=0 
      JVAL=JVAL+1 
      JZER=JZER+1 
32508 CONTINUE 
756   Format(' ICN LUDI:',10I6) 
32511 CONTINUE 
32694 CONTINUE 
32693 CONTINUE 
      IF(ILAST.NE.NN)IACTIV=IPTR(ILAST+1) 
32744 CONTINUE 
32751 CONTINUE 
      IF(.NOT.(IRANK.NE.NN))GOTO32505 
      DO32502I=1,NN 
      IF(.NOT.(LASTC(I).LE.0))GOTO32499 
      ISING=-LASTC(I) 
      IQ(ISING)=-IQ(ISING) 
      LASTC(I)=ISING 
32499 CONTINUE 
32502 CONTINUE 
32505 CONTINUE 
      ISTART=IDISP(1) 
      IEND=IBEG-1 
      DO32496JJ=ISTART,IEND 
      JOLD=ICN(JJ) 
      ICN(JJ)=LASTC(JOLD) 
32496 CONTINUE 
      DO32493II=1,NN 
      I=LASTR(II) 
      NEXTR(I)=LENR(II) 
      NEXTC(I)=LENRL(II) 
32493 CONTINUE 
      DO32490I=1,NN 
      LENRL(I)=NEXTC(I) 
      LENR(I)=NEXTR(I) 
32490 CONTINUE 
      DO32487II=1,NN 
      I=LASTR(II) 
      J=LASTC(II) 
      NEXTR(I)=IABS(IP(II)) 
      NEXTC(J)=IABS(IQ(II)) 
32487 CONTINUE 
      DO32484I=1,NN 
      IF(IP(I).LT.0)NEXTR(I)=-NEXTR(I) 
      IP(I)=NEXTR(I) 
      IF(IQ(I).LT.0)NEXTC(I)=-NEXTC(I) 
      IQ(I)=NEXTC(I) 
32484 CONTINUE 
      IP(NN)=IABS(IP(NN)) 
      IDISP(2)=IEND 
      IF(IFAIL.NE.5)GOTO2060 
      IF(MOD(ISAVE/10,10).NE.0)WRITE(NERR,99991)MINICN 
      GOTO2060 
1980  continue 
      IF(.NOT.(IFAIL.NE.5))GOTO32481 
      IFAIL=3 
      GOTO32480 
32481 CONTINUE 
      IFAIL=6 
32480 CONTINUE 
2020  CONTINUE 
      IF(.NOT.(MOD(ISAVE/10,10).NE.0))GOTO32478 
      GOTO(32475,32474,32473,32472,32471,32470),IFAIL 
      GOTO32469 
32475 CONTINUE 
      WRITE(NERR,99999) 
      GOTO32469 
32474 CONTINUE 
      WRITE(NERR,99998) 
      GOTO32469 
32473 CONTINUE 
      WRITE(NERR,99996) 
      GOTO32469 
32472 CONTINUE 
      WRITE(NERR,99995) 
      GOTO32469 
32471 CONTINUE 
      WRITE(NERR,99994) 
      GOTO32469 
32470 CONTINUE 
      WRITE(NERR,99993) 
32469 CONTINUE 
      IF(.NOT.(IFAIL.GT.2))GOTO32468 
      PIVOT=PIVOT-ISTART+1 
      WRITE(NERR,99992)PIVOT,NBLOCK,ISTART,ILAST 
      IF(PIVOT.EQ.0)WRITE(NERR,99990)MINIRN 
32468 CONTINUE 
32478 CONTINUE 
      IDISP(2)=IACTIV 
2060  IDISP(3)=IRNCP 
      IDISP(4)=ICNCP 
      IDISP(5)=IRANK 
      IDISP(6)=MINIRN 
      IDISP(7)=MINICN 
      RETURN 
99999 FORMAT(50H0MATRIX IS STRUCTURALLY SINGULAR - DECOMPOSITION A,6HBOR 
     *TED) 
99998 FORMAT(50H0MATRIX IS NUMERICALLY SINGULAR - DECOMPOSITION AB,5HORT 
     *ED) 
99997 FORMAT(48H0LU DECOMPOSITION DESTROYED TO CREATE MORE SPACE) 
99996 FORMAT(17H0LIRN TOO SMALL -) 
99995 FORMAT(22H0LICN MUCH TOO SMALL -) 
99994 FORMAT(17H0LICN TOO SMALL -) 
99993 FORMAT(26H0LICN AND LIRN TOO SMALL -) 
99992 FORMAT(33H  DECOMPOSITION ABORTED AT STAGE ,I5,9H IN BLOCK,1H ,I5/ 
     *17H  WITH FIRST ROW ,I5,14H AND LAST ROW ,I5) 
99991 FORMAT(17H0LICN TOO SMALL -/29H  FOR SUCCESSFUL DECOMPOSITIO,23HN  
     *SET LICN TO AT LEAST ,I8) 
99990 FORMAT(35H  TO CONTINUE SET LIRN TO AT LEAST ,I8) 
      END 

      SUBROUTINEM28BYS(N,NZ,A,LICN,IVECT,JVECT,ICN,IKEEP,IW,W,GROW,EPS,R 
     *MIN,ABORT,IDISP,IFAIL) 
      DOUBLEPRECISIONEPS,RMIN 
      INTEGERIFAIL,LICN,N,NZ 
      LOGICALABORT,GROW 
      DOUBLEPRECISIONA(LICN),W(N) 
      INTEGERICN(LICN),IDISP(2),IKEEP(N,5),IVECT(NZ),IW(N,5),JVECT(NZ) 
      CHARACTER*6SRNAME 
      PARAMETER(SRNAME='M28BYS') 
      DOUBLEPRECISIONWMAX 
      INTEGERI1,IEND,IR,ISAVE,NERR 
      CHARACTER*1P01REC(1) 
      CHARACTER*65REC(2) 
      INTEGERP01ABF 
      EXTERNALP01ABF 
c      EXTERNALF01BRQ,M30BYS,F01DYS,X04AAF,X04BAF !!!!!!!!!!!!!!!!!!!!
      EXTERNALF01BRQ,M30BYS,       X04AAF,X04BAF 
      INTRINSICMOD 
      ISAVE=IFAIL 
      CALLX04AAF(0,NERR) 
      IF(N.GT.0)GOTO20 
      IFAIL=1 
      GOTO120 
20    IF(NZ.GT.0)GOTO40 
      IFAIL=2 
      GOTO120 
40    IF(LICN.GE.NZ)GOTO60 
      IFAIL=3 
      GOTO120 
60    CALLM28DYS(N,A,LICN,IVECT,JVECT,NZ,ICN,IKEEP,IKEEP(1,4),IKEEP(1,5) 
     *,IKEEP(1,2),IKEEP(1,3),IW(1,3),IW,W(1),ABORT,IDISP,IFAIL) 
      WMAX=W(1) 
      IF(IFAIL.NE.0)GOTO120 
      IFAIL=ISAVE 
      CALLM30BYS(N,ICN,A,LICN,IKEEP,IKEEP(1,4),IDISP,IKEEP(1,2),IKEEP(1, 
     *3),W,IW,EPS,RMIN,IFAIL) 
      IF(IFAIL.EQ.0)GOTO100 
      IF(IFAIL.GT.0)GOTO80 
      IR=-IFAIL 
      IFAIL=6 
      GOTO120 
80    IR=IFAIL 
      IFAIL=7 
100   I1=IDISP(1) 
      IEND=LICN-I1+1 
      IF(GROW)CALLF01BRQ(N,ICN,A(I1),IEND,IKEEP,IKEEP(1,4),W) 
      IF(GROW)W(1)=W(1)+WMAX 
      IF(IFAIL.EQ.0)GOTO160 
120   CONTINUE 
      IF(MOD(ISAVE/10,10).EQ.0)GOTO140 
      CALLX04AAF(0,NERR) 
      GOTO(32757,32756,32755,32754,32753,32752,32751,32750),IFAIL 
      GOTO32749 
32757 CONTINUE 
      WRITE(REC,FMT=99999)N 
      GOTO32749 
32756 CONTINUE 
      WRITE(REC,FMT=99998)NZ 
      GOTO32749 
32755 CONTINUE 
      WRITE(REC,FMT=99997)LICN,NZ 
      GOTO32749 
32754 CONTINUE 
      GOTO32749 
32753 CONTINUE 
      GOTO32749 
32752 CONTINUE 
      WRITE(REC,FMT=99996)IR 
      IW(1,1)=IR 
      GOTO32749 
32751 CONTINUE 
      WRITE(REC,FMT=99995)IR 
      GOTO32749 
32750 CONTINUE 
      WRITE(REC,99994) 
32749 CONTINUE 
      IF((IFAIL.GE.1.AND.IFAIL.LE.3).OR.(IFAIL.GE.6.AND.IFAIL.LE.8))THEN 
      CALLX04BAF(NERR,REC(1)) 
      CALLX04BAF(NERR,REC(2)) 
      ENDIF 
      IF(IFAIL.EQ.7)GOTO160 
140   IFAIL=P01ABF(ISAVE,IFAIL,SRNAME,0,P01REC) 
160   RETURN 
99999 FORMAT(/' ON ENTRY N .LE. 0 , N =',I10) 
99998 FORMAT(/' ON ENTRY NZ .LE. 0 , NZ =',I10) 
99997 FORMAT(/' ON ENTRY LICN .LT. NZ , LICN = ',I8,'  NZ = ',I8) 
99996 FORMAT(/' NUMERICAL SINGULARITY IN ROW ',I4,' - DECOMPOSITION A',' 
     *BORTED') 
99995 FORMAT(/' SUBTHRESHOLD PIVOT IN ROW ',I4,' - DECOMPOSITION COMP',' 
     *LETED') 
99994 FORMAT(/' DUPLICATE ELEMENTS FOUND ON INPUT - SEE ADVISORY MESS',' 
     *AGES') 
      END 
      SUBROUTINEM30BYS(N,ICN,A,LICN,LENR,LENRL,IDISP,IP,IQ,W,IW,EPS,RMIN 
     *,IFAIL) 
      DOUBLEPRECISIONEPS,RMIN 
      INTEGERIFAIL,LICN,N 
      DOUBLEPRECISIONA(LICN),W(N) 
      INTEGERICN(LICN),IDISP(2),IP(N),IQ(N),IW(N,2),LENR(N),LENRL(N) 
      DOUBLEPRECISIONAU,ONE,ROWMAX,ZERO 
      INTEGERI,IFIN,ILEND,IPIVJ,ISING,ISTART,J,JAY,JAYJAY,JFIN,JJ,PIVPOS 
     *,KIW 
      LOGICALSTAB 
      DATAZERO/0.0D0/,ONE/1.0D0/ 
      STAB=EPS.LE.ONE 
      RMIN=EPS 
      ISING=0 
      IFAIL=0 
      IF(N.EQ.1)GOTO320 
      KIW=0 
      DO20I=1,N 
      W(I)=ZERO 
20    CONTINUE 
      IW(1,1)=IDISP(1) 
      DO40I=2,N 
      IW(I,1)=IW(I-1,1)+LENR(I-1) 
40    CONTINUE 
      DO280I=1,N 
      ISTART=IW(I,1) 
      IFIN=ISTART+LENR(I)-1 
      ILEND=ISTART+LENRL(I)-1 
      IF(ISTART.GT.ILEND)GOTO140 
      DO60JJ=ISTART,IFIN 
      J=ICN(JJ) 
      W(J)=A(JJ) 
60    CONTINUE 
      DO100JJ=ISTART,ILEND 
      J=ICN(JJ) 
      IPIVJ=IW(J,1)+LENRL(J) 
      AU=-W(J)/A(IPIVJ) 
      IF(.NOT.(W(J).EQ.ZERO))GOTO32748 
      KIW=KIW+1 
      IW(KIW,2)=J 
32748 CONTINUE 
      W(J)=AU 
      IPIVJ=IPIVJ+1 
      JFIN=IW(J,1)+LENR(J)-1 
      IF(IPIVJ.GT.JFIN)GOTO100 
      DO80JAYJAY=IPIVJ,JFIN 
      JAY=ICN(JAYJAY) 
      IF(.NOT.(W(JAY).EQ.ZERO))GOTO32745 
      KIW=KIW+1 
      IW(KIW,2)=JAY 
32745 CONTINUE 
      W(JAY)=W(JAY)+AU*A(JAYJAY) 
80    CONTINUE 
100   CONTINUE 
      DO120JJ=ISTART,IFIN 
      J=ICN(JJ) 
      A(JJ)=W(J) 
120   CONTINUE 
      DO32742JAY=1,KIW 
      J=IW(JAY,2) 
      W(J)=ZERO 
32742 CONTINUE 
      KIW=0 
140   PIVPOS=ILEND+1 
      IF(IQ(I).GT.0)GOTO240 
      IF(ISING.EQ.0)ISING=I 
      IF(PIVPOS.GT.IFIN)GOTO160 
      IF(A(PIVPOS).NE.ZERO)GOTO300 
160   IF(ISTART.GT.IFIN)GOTO200 
      DO180JJ=ISTART,IFIN 
      IF(ICN(JJ).LT.ISING)GOTO180 
      IF(A(JJ).NE.ZERO)GOTO300 
180   CONTINUE 
200   IF(PIVPOS.LE.IFIN)A(PIVPOS)=ONE 
      IF(IP(I).GT.0.AND.I.NE.N)GOTO280 
      DO220J=ISING,I 
      IF((LENR(J)-LENRL(J)).EQ.0)GOTO220 
      JJ=IW(J,1)+LENRL(J) 
      A(JJ)=ZERO 
220   CONTINUE 
      ISING=0 
      GOTO280 
240   IF(PIVPOS.GT.IFIN)GOTO300 
      IF(A(PIVPOS).EQ.ZERO)GOTO300 
      IF(.NOT.STAB)GOTO280 
      ROWMAX=ZERO 
      DO260JJ=PIVPOS,IFIN 
      ROWMAX=DMAX1(ROWMAX,DABS(A(JJ))) 
260   CONTINUE 
      IF(DABS(A(PIVPOS))/ROWMAX.GE.RMIN)GOTO280 
      IFAIL=I 
      RMIN=DABS(A(PIVPOS))/ROWMAX 
280   CONTINUE 
      GOTO320 
300   IFAIL=-I 
320   RETURN 
      END 
      SUBROUTINEM28DYS(N,A,LICN,IVECT,JVECT,NZ,ICN,LENR,LENRL,LENOFF,IP, 
     *IQ,IW1,IW,W1,ABORT,IDISP,IFAIL) 
      DOUBLEPRECISIONW1 
      INTEGERIFAIL,LICN,N,NZ 
      LOGICALABORT 
      DOUBLEPRECISIONA(LICN) 
      INTEGERICN(LICN),IDISP(2),IP(N),IQ(N),IVECT(NZ),IW1(N,3),IW(N,2),J 
     *VECT(NZ),LENOFF(N),LENR(N),LENRL(N) 
      DOUBLEPRECISIONAA,ZERO 
      INTEGERI,IBLOCK,IDISP2,IDUMMY,II,INEW,IOLD,ISAVE,J1,J2,JCOMP,JDUMM 
     *Y,JJ,JNEW,JOLD,MIDPT,NADV,NERR 
      LOGICALBLOCKL 
      DATAZERO/0.0D0/ 
      ISAVE=IFAIL 
      IFAIL=0 
      CALLX04AAF(0,NERR) 
      CALLX04ABF(0,NADV) 
      BLOCKL=LENOFF(1).GE.0 
      IBLOCK=1 
      IW(1,1)=1 
      IW(1,2)=IDISP(1) 
      DO20I=1,N 
      IW1(I,3)=IBLOCK 
      IF(IP(I).LT.0)IBLOCK=IBLOCK+1 
      II=IABS(IP(I)) 
      IW1(II,1)=I 
      JJ=IQ(I) 
      JJ=IABS(JJ) 
      IW1(JJ,2)=I 
      IF(I.EQ.1)GOTO20 
      IF(BLOCKL)IW(I,1)=IW(I-1,1)+LENOFF(I-1) 
      IW(I,2)=IW(I-1,2)+LENR(I-1) 
20    CONTINUE 
      IDISP2=IDISP(2) 
      DO340I=1,NZ 
      IF(I.GT.IDISP2)GOTO40 
      IF(ICN(I).LT.0)GOTO340 
40    IOLD=IVECT(I) 
      JOLD=JVECT(I) 
      AA=A(I) 
      DO280IDUMMY=1,NZ 
      IF(IOLD.LE.N.AND.IOLD.GT.0.AND.JOLD.LE.N.AND.JOLD.GT.0)GOTO60 
      IFAIL=4 
      IF(MOD(ISAVE/10,10).NE.0)WRITE(NERR,99999)I,A(I),IOLD,JOLD 
      GOTO360 
60    INEW=IW1(IOLD,1) 
      JNEW=IW1(JOLD,2) 
      IF(IW1(INEW,3)-IW1(JNEW,3))80,120,100 
80    IFAIL=5 
      IF(MOD(ISAVE/10,10).NE.0)WRITE(NERR,99998)IOLD,JOLD 
      GOTO360 
100   J1=IW(INEW,1) 
      J2=J1+LENOFF(INEW)-1 
      GOTO220 
120   J1=IW(INEW,2) 
      IF(INEW.GT.JNEW)GOTO140 
      J2=J1+LENR(INEW)-1 
      J1=J1+LENRL(INEW) 
      GOTO220 
140   J2=J1+LENRL(INEW) 
      DO200JDUMMY=1,N 
      MIDPT=(J1+J2)/2 
      JCOMP=IABS(ICN(MIDPT)) 
      IF(JNEW-JCOMP)160,260,180 
160   J2=MIDPT 
      GOTO200 
180   J1=MIDPT 
200   CONTINUE 
      IFAIL=5 
      IF(MOD(ISAVE/10,10).NE.0)WRITE(NERR,99997)IOLD,JOLD 
      GOTO360 
220   DO240MIDPT=J1,J2 
      IF(IABS(ICN(MIDPT)).EQ.JNEW)GOTO260 
240   CONTINUE 
      IFAIL=5 
      IF(MOD(ISAVE/10,10).NE.0)WRITE(NERR,99997)IOLD,JOLD 
      GOTO360 
260   IF(ICN(MIDPT).LT.0)GOTO320 
      IF(MIDPT.GT.NZ.OR.MIDPT.LE.I)GOTO300 
      W1=A(MIDPT) 
      A(MIDPT)=AA 
      AA=W1 
      IOLD=IVECT(MIDPT) 
      JOLD=JVECT(MIDPT) 
      ICN(MIDPT)=-ICN(MIDPT) 
280   CONTINUE 
300   A(MIDPT)=AA 
      ICN(MIDPT)=-ICN(MIDPT) 
      GOTO340 
320   A(MIDPT)=A(MIDPT)+AA 
      IF(MOD(ISAVE/100,100).NE.0)WRITE(NADV,99996)IOLD,JOLD,AA 
      IF(ABORT)IFAIL=8 
340   CONTINUE 
360   W1=ZERO 
      DO400I=1,IDISP2 
      IF(ICN(I).LT.0)GOTO380 
      A(I)=ZERO 
      GOTO400 
380   ICN(I)=-ICN(I) 
      W1=DMAX1(W1,DABS(A(I))) 
400   CONTINUE 
      RETURN 
99999 FORMAT(49H0ON ENTRY IRN(I) OR ICN(I) IS OUT OF RANGE - I = ,I8/9H  
     * A(I) = ,1PD14.6,11H  IRN(I) = ,I8,11H  ICN(I) = ,I8) 
99998 FORMAT(19H0NON-ZERO ELEMENT (,I6,1H,,I6,13H) IN ZERO OFF,15H-DIAGO 
     *NAL BLOCK) 
99997 FORMAT(19H0NON-ZERO ELEMENT (,I6,1H,,I6,13H) WAS NOT IN ,11HL/U PA 
     *TTERN) 
99996 FORMAT(45H0M28DYS FOUND DUPLICATE ELEMENT WITH INDICES ,I6,2H, ,I6 
     */10H  VALUE = ,1PD14.6) 
      END 

      SUBROUTINEM28CYS(N,A,LICN,ICN,IKEEP,RHS,W,MTYPE,IDISP,RESID) 
      DOUBLEPRECISIONRESID 
      INTEGERLICN,MTYPE,N 
      DOUBLEPRECISIONA(LICN),RHS(N),W(N) 
      INTEGERICN(LICN),IDISP(2),IKEEP(N,5) 
      CALLM30CYS(N,ICN,A,LICN,IKEEP,IKEEP(1,4),IKEEP(1,5),IDISP,IKEEP(1, 
     *2),IKEEP(1,3),RHS,W,MTYPE,RESID) 
      RETURN 
      END 
      SUBROUTINEM30CYS(N,ICN,A,LICN,LENR,LENRL,LENOFF,IDISP,IP,IQ,X,W,MT 
     *YPE,RESID) 
      DOUBLEPRECISIONRESID 
      INTEGERLICN,MTYPE,N 
      DOUBLEPRECISIONA(LICN),W(N),X(N) 
      INTEGERICN(LICN),IDISP(2),IP(N),IQ(N),LENOFF(N),LENR(N),LENRL(N) 
      DOUBLEPRECISIONWI,WII,ZERO 
      INTEGERI,IB,IBACK,IBLEND,IBLOCK,IEND,IFIRST,II,III,ILAST,J1,J2,J3, 
     *J,JJ,JPIV,JPIVP1,K,LJ1,LJ2,LT,LTEND,NUMBLK 
      LOGICALNEG,NOBLOC 
      DATAZERO/0.0D0/ 
      RESID=ZERO 
      NOBLOC=LENOFF(1).LT.0 
      IF(MTYPE.NE.1)GOTO280 
      NEG=.FALSE. 
      IP(N)=-IP(N) 
      DO20II=1,N 
      I=IP(II) 
      I=IABS(I) 
      W(II)=X(I) 
20    CONTINUE 
      LT=1 
      IFIRST=1 
      IBLOCK=IDISP(1) 
      DO240I=1,N 
      WI=W(I) 
      IF(NOBLOC)GOTO60 
      IF(LENOFF(I).EQ.0)GOTO60 
      LTEND=LT+LENOFF(I)-1 
      DO40JJ=LT,LTEND 
      J=ICN(JJ) 
      WI=WI-A(JJ)*W(J) 
40    CONTINUE 
      LT=LTEND+1 
60    IF(IP(I).LT.0)NEG=.TRUE. 
      IF(LENRL(I).EQ.0)GOTO100 
      IEND=IBLOCK+LENRL(I)-1 
      DO80JJ=IBLOCK,IEND 
      J=ICN(JJ) 
      WI=WI+A(JJ)*W(J) 
80    CONTINUE 
100   IBLOCK=IBLOCK+LENR(I) 
      W(I)=WI 
      IF(.NOT.NEG)GOTO240 
      J1=IBLOCK 
      IB=I 
      IF(IQ(I).GT.0)GOTO140 
      DO120III=IFIRST,I 
      IB=I-III+IFIRST 
      IF(IQ(IB).GT.0)GOTO140 
      J1=J1-LENR(IB) 
      RESID=DMAX1(RESID,DABS(W(IB))) 
      W(IB)=ZERO 
120   CONTINUE 
      GOTO220 
140   DO200III=IFIRST,IB 
      II=IB-III+IFIRST 
      J2=J1-1 
      J1=J1-LENR(II) 
      JPIV=J1+LENRL(II) 
      JPIVP1=JPIV+1 
      IF(J2.LT.JPIVP1)GOTO180 
      WII=W(II) 
      DO160JJ=JPIVP1,J2 
      J=ICN(JJ) 
      WII=WII-A(JJ)*W(J) 
160   CONTINUE 
      W(II)=WII 
180   W(II)=W(II)/A(JPIV) 
200   CONTINUE 
220   IFIRST=I+1 
      NEG=.FALSE. 
240   CONTINUE 
      DO260II=1,N 
      I=IQ(II) 
      I=IABS(I) 
      X(I)=W(II) 
260   CONTINUE 
      IP(N)=-IP(N) 
      GOTO640 
280   DO300II=1,N 
      I=IQ(II) 
      I=IABS(I) 
      W(II)=X(I) 
300   CONTINUE 
      LJ1=IDISP(1) 
      IBLOCK=IDISP(2)+1 
      ILAST=N 
      IBLEND=IBLOCK 
      DO580NUMBLK=1,N 
      IF(ILAST.EQ.0)GOTO600 
      IBLOCK=IBLOCK-LENR(ILAST) 
      DO320K=1,N 
      II=ILAST-K 
      IF(II.EQ.0)GOTO340 
      IF(IP(II).LT.0)GOTO340 
      IBLOCK=IBLOCK-LENR(II) 
320   CONTINUE 
340   IFIRST=II+1 
      J1=IBLOCK 
      DO420I=IFIRST,ILAST 
      IF(W(I).EQ.ZERO)GOTO400 
      IF(IQ(I).LT.0)GOTO440 
      J2=J1+LENRL(I) 
      WI=W(I)/A(J2) 
      IF(LENR(I)-LENRL(I).EQ.1)GOTO380 
      J2=J2+1 
      J3=J1+LENR(I)-1 
      DO360JJ=J2,J3 
      J=ICN(JJ) 
      W(J)=W(J)-A(JJ)*WI 
360   CONTINUE 
380   W(I)=WI 
400   J1=J1+LENR(I) 
420   CONTINUE 
      GOTO480 
440   DO460II=I,ILAST 
      RESID=DMAX1(RESID,DABS(W(II))) 
      W(II)=ZERO 
460   CONTINUE 
480   J1=IBLEND 
      DO560IBACK=IFIRST,ILAST 
      I=ILAST-IBACK+IFIRST 
      J1=J1-LENR(I) 
      IF(LENRL(I).EQ.0)GOTO520 
      J2=J1+LENRL(I)-1 
      DO500JJ=J1,J2 
      J=ICN(JJ) 
      W(J)=W(J)+A(JJ)*W(I) 
500   CONTINUE 
520   IF(NOBLOC)GOTO560 
      IF(LENOFF(I).EQ.0)GOTO560 
      LJ2=LJ1-1 
      LJ1=LJ1-LENOFF(I) 
      DO540JJ=LJ1,LJ2 
      J=ICN(JJ) 
      W(J)=W(J)-A(JJ)*W(I) 
540   CONTINUE 
560   CONTINUE 
      IBLEND=J1 
      ILAST=IFIRST-1 
580   CONTINUE 
600   DO620II=1,N 
      I=IP(II) 
      I=IABS(I) 
      X(I)=W(II) 
620   CONTINUE 
640   RETURN 
      END 

      SUBROUTINEM28CYN(N,A,LICN,ICN,IKEEP,RHS,W,IW,MTYPE,IDISP,RESID,Ifa 
     *il) 
      DOUBLEPRECISIONRESID 
      INTEGERLICN,MTYPE,N 
      DOUBLEPRECISIONA(LICN),RHS(N),W(N) 
      INTEGERICN(LICN),IDISP(2),IKEEP(N,5) 
      CALLM30CYN(N,ICN,A,LICN,IKEEP,IKEEP(1,4),IKEEP(1,5),IDISP,IKEEP(1, 
     *2),IKEEP(1,3),RHS,W,MTYPE,RESID,Ifail) 
      RETURN 
      END 
      SUBROUTINEM30CYN(N,ICN,A,LICN,LENR,LENRL,LENOFF,IDISP,IP,IQ,X,W,MT 
     *YPE,RESID,Ifail) 
      DOUBLEPRECISIONRESID 
      INTEGERLICN,MTYPE,N 
      DOUBLEPRECISIONA(LICN),W(N),X(N) 
      INTEGERICN(LICN),IDISP(2),IP(N),IQ(N),LENOFF(N),LENR(N),LENRL(N) 
      DOUBLEPRECISIONWI,WII,ZERO 
      INTEGERI,IB,IBACK,IBLEND,IBLOCK,IEND,IFIRST,II,III,ILAST,J1,J2,J3, 
     *J,JJ,JPIV,JPIVP1,K,LJ1,LJ2,LT,LTEND,NUMBLK 
      LOGICALNEG,NOBLOC 
      DATAZERO/0.0D0/ 
      Ifail=0 
      RESID=ZERO 
      NOBLOC=LENOFF(1).LT.0 
      IF(MTYPE.NE.1)GOTO280 
      NEG=.FALSE. 
      IP(N)=-IP(N) 
      DO20II=1,N 
      I=IP(II) 
      I=IABS(I) 
      W(II)=X(I) 
20    CONTINUE 
      LT=1 
      IFIRST=1 
      IBLOCK=IDISP(1) 
      DO240I=1,N 
      WI=W(I) 
      IF(.NOT.NOBLOC.AND.LENOFF(I).NE.0)Then 
      LTEND=LT+LENOFF(I)-1 
      DO40JJ=LT,LTEND 
      J=ICN(JJ) 
      If(abs(A(JJ)).LT.1.D0)Then 
      WI=WI-A(JJ)*W(J) 
      elseIf(abs(W(J)).LT.1.D+70/abs(A(JJ)))Then 
      WI=WI-A(JJ)*W(J) 
      else 
      Ifail=13 
      Return 
      endif 
40    CONTINUE 
      LT=LTEND+1 
      Endif 
      IF(IP(I).LT.0)NEG=.TRUE. 
      IF(LENRL(I).NE.0)Then 
      IEND=IBLOCK+LENRL(I)-1 
      DO80JJ=IBLOCK,IEND 
      J=ICN(JJ) 
      If(abs(A(JJ)).LT.1.D0)Then 
      WI=WI+A(JJ)*W(J) 
      elseIf(abs(W(J)).LT.1.D+70/abs(A(JJ)))Then 
      WI=WI+A(JJ)*W(J) 
      else 
      Ifail=13 
      Return 
      endif 
80    CONTINUE 
      Endif 
      IBLOCK=IBLOCK+LENR(I) 
      W(I)=WI 
      IF(.NOT.NEG)GOTO240 
      J1=IBLOCK 
      IB=I 
      IF(IQ(I).GT.0)GOTO140 
      DO120III=IFIRST,I 
      IB=I-III+IFIRST 
      IF(IQ(IB).GT.0)GOTO140 
      J1=J1-LENR(IB) 
      RESID=DMAX1(RESID,DABS(W(IB))) 
      W(IB)=ZERO 
120   CONTINUE 
      GOTO220 
140   DO200II=IB,IFIRST,-1 
      J2=J1-1 
      J1=J1-LENR(II) 
      JPIV=J1+LENRL(II) 
      JPIVP1=JPIV+1 
      IF(J2.GE.JPIVP1)Then 
      WII=W(II) 
      DO160JJ=JPIVP1,J2 
      J=ICN(JJ) 
      WII=WII-A(JJ)*W(J) 
160   CONTINUE 
      W(II)=WII 
      Endif 
      If(1.D+35.LT.abs(A(JPIV)))Then 
      W(II)=W(II)/A(JPIV) 
      elseIf(abs(W(II)).LT.1.D+40*abs(A(JPIV)))Then 
      W(II)=W(II)/A(JPIV) 
      else 
      Ifail=13 
      Return 
      endif 
200   CONTINUE 
220   IFIRST=I+1 
      NEG=.FALSE. 
240   CONTINUE 
      DO260II=1,N 
      I=IQ(II) 
      I=IABS(I) 
      X(I)=W(II) 
260   CONTINUE 
      IP(N)=-IP(N) 
      GOTO640 
280   DO300II=1,N 
      I=IQ(II) 
      I=IABS(I) 
      W(II)=X(I) 
300   CONTINUE 
      LJ1=IDISP(1) 
      IBLOCK=IDISP(2)+1 
      ILAST=N 
      IBLEND=IBLOCK 
      DO580NUMBLK=1,N 
      IF(ILAST.EQ.0)GOTO600 
      IBLOCK=IBLOCK-LENR(ILAST) 
      DO320K=1,N 
      II=ILAST-K 
      IF(II.EQ.0)GOTO340 
      IF(IP(II).LT.0)GOTO340 
      IBLOCK=IBLOCK-LENR(II) 
320   CONTINUE 
340   IFIRST=II+1 
      J1=IBLOCK 
      DO420I=IFIRST,ILAST 
      IF(W(I).EQ.ZERO)GOTO400 
      IF(IQ(I).LT.0)GOTO440 
      J2=J1+LENRL(I) 
      WI=W(I)/A(J2) 
      IF(LENR(I)-LENRL(I).EQ.1)GOTO380 
      J2=J2+1 
      J3=J1+LENR(I)-1 
      DO360JJ=J2,J3 
      J=ICN(JJ) 
      W(J)=W(J)-A(JJ)*WI 
360   CONTINUE 
380   W(I)=WI 
400   J1=J1+LENR(I) 
420   CONTINUE 
      GOTO480 
440   DO460II=I,ILAST 
      RESID=DMAX1(RESID,DABS(W(II))) 
      W(II)=ZERO 
460   CONTINUE 
480   J1=IBLEND 
      DO560IBACK=IFIRST,ILAST 
      I=ILAST-IBACK+IFIRST 
      J1=J1-LENR(I) 
      IF(LENRL(I).EQ.0)GOTO520 
      J2=J1+LENRL(I)-1 
      DO500JJ=J1,J2 
      J=ICN(JJ) 
      W(J)=W(J)+A(JJ)*W(I) 
500   CONTINUE 
520   IF(NOBLOC)GOTO560 
      IF(LENOFF(I).EQ.0)GOTO560 
      LJ2=LJ1-1 
      LJ1=LJ1-LENOFF(I) 
      DO540JJ=LJ1,LJ2 
      J=ICN(JJ) 
      W(J)=W(J)-A(JJ)*W(I) 
540   CONTINUE 
560   CONTINUE 
      IBLEND=J1 
      ILAST=IFIRST-1 
580   CONTINUE 
600   DO620II=1,N 
      I=IP(II) 
      I=IABS(I) 
      X(I)=W(II) 
620   CONTINUE 
640   RETURN 
      END 
